Pattern Dynamics of Vortex Ripples in Sand: 
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Vortex ripples in sand are studied experimentally in a one-dimensional setup with periodic 
boundary conditions. The nonlinear evolution, far from the onset of instability, is analyzed in 
the framework of a simple model developed for homogeneous patterns. The interaction function 
describing the mass transport between neighboring ripples is extracted from experimental runs 
using a recently proposed method for data analysis, and the predictions of the model are compared 
to the experiment. An analytic explanation of the wavelength selection mechanism in the model is 
provided, and the width of the stable band of ripples is measured. 



Ever since the establishment of a conceptual frame- 
work for pattern formation [[l] , the description of patterns 
formed in sand by the flow of wind or water has posed a 
challenge to the community . Despite diverse efforts 
including coupled map models for ripples and dunes in air 
stochastic models for ripples in air || or continuum 
equations based on the symmetries of the problem |Q, 
theoretical understanding has remained sparse. For ex- 
ample, all models display a coarsening of the ripple/dune 
pattern, but the coarsening does not terminate at a final 
selected wavelength, as is frequently observed in nature. 
Furthermore the models are heuristic, and it is not possi- 
ble to make a quantitative comparison with experiments. 

Here, we study vortex ripples which are created 
by an oscillatory water flow, such as that generated near 
the sand bed by a surface wave. Vortex ripples have at- 
tracted attention as an example of a non-linear pattern 
forming system with a strongly sub-critical first bifurca- 
tion p,|7|,|9|j which cannot be described by conventional 
methods like amplitude equations |j) . As most other sand 
patterns, they display coarsening and saturation at a fi- 
nite wavelength. The approach pursued in this Letter 
combines a simple model for the fully developed pattern 
with a sophisticated data analysis which allows to extract 
the key model ingredient - the interaction function /(A) 
- directly from the experimental runs. In this way the 
validity of the model can be tested, and additional fea- 
tures required for the description, can be identified. As 
far as we know, this is the first quantitative comparison 
between theory and experiment for a sand pattern. The 
basic ideas are general, and we expect that the theoret- 
ical formalism combined with the experimental analysis 
can be used for related sand patters (e.g. the one studied 
in j^]), or other strongly nonlinear systems. 

One dimensional vortex ripples can be created in an 
annular channel |^|,|l0|, ensuring the pattern to be sub- 
ject to well-defined, periodic boundary conditions. Freely 



grown ripples are created from a flat bed by a coarsen- 
ing process, which eventually saturates at an equilibrium 
state, where the ripple length A is almost independent of 
the frequency v of the driving and proportional to the 
amplitude of the oscillation of the plate a |^,[ll| • Here we 
are mostly concerned with the stability and evolution of 
the ripple patterns themselves, and not with the insta- 
bility of the flat bed, which has been discussed elsewhere 
p2|,p|. By creating a homogeneous ripple pattern, where 
all ripples have the same length, and changing amplitude 
and frequency, the stability of the pattern is probed. In 
this way it was found that there is a stable band of rip- 
ples, X m in < A < X m ax, for a given set of driving conditions 

0- 

Recently, a simple model was proposed, which de- 
scribes both the coarsening and saturation of ripples, and 
reproduces the existence of a stable band ||. In its sim- 
plest form, the change of the length A.,- of the ripple j is 
a function of Xj itself and the lengths of the neighboring 
ripples, 

A,- = -/(A i _i) + 2/(A i )-/(A i+ i) 1 (1) 

where the interaction function /(A) describes the trans- 
fer of mass between neighbouring ripples (Fig. [j]). Argu- 
ments based on the properties of the separation vortex 
forming behind a ripple show that /(A) should be convex, 
its approximate shape has been found from a numerical 
simulation of the fluid and sand flow over homogeneous 
ripples M. The purpose of this Letter is a test of the 
model (El) through comparison with experimental data. 

Experimental setup. All mm wide, 15 cm high annu- 
lar plexiglass channel of diameter 48.6 cm is filled with 
water, A, and glass beads, B, with a diameter of 250 
±50 fj,m (Fig. ||). In the middle is a conical mirror, C, 
which is filmed from above by a CCD camera with a res- 
olution of 640x480 pixels. The channel is attached to a 
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motor, D, by a 1.5 m long arm E, which oscillates the 
channel in an almost sinusoidal fashion. 

We create initial conditions with small ripples by oscil- 
lating at a small amplitude. The experiment starts when 
amplitude and frequency are changed to their desired val- 
ues. Here, we focus on the condition a = 6 cm, v = 0.6 
Hz. We want to emphasize that there is nothing special 
about this condition, the same qualitative results are ob- 
tained for other parameters. Growing the ripples from 
initial conditions with homogeneous ripples smaller than 
Amm results in a final number of 19-21 ripples with a 
length (averaged over 16 realizations) of A/a=1.25. The 
uncertainty in the length due to periodicity is around 
5 %. The images are transformed back to cartesian coor- 
dinates, and the lengths Xj are found by fitting them to 
triangles with fixed slope (Fig. ||a). From series of con- 
secutive Xj, the temporal change, Xj, is established, and 
a space-time plot of the ripple evolution is constructed. 
Fig. [|b shows the evolution starting from a small initial 
wavelength, which leads to the annihilation of ripples. 
Around each annihilation the ripples have not been anal- 
ysed, as it becomes inaccurate to fit the triangles during 
these events. 

Measuring the stable band. By running the experi- 
ment with different values of a until a homogeneous pat- 
tern occurs, we create initial conditions with different 
no. of ripples. From this initial pattern, the experi- 
ment is started with the above parameters, and run for 
10000 periods. We have created initial conditions with 
N e [17; 24] and observed that for N > 22 ripples are 
annihilated, while for N < 18 one or more new ripples 
are created. We therefore conclude that there exists a 
stable band N G [17.5; 22.5] (for a = 6 cm, v = 0.6 
Hz), which corresponds to X m i n /a = 1.13 ± 0.03 and 
Xmax/o, — 1.45 ± 0.05. Note that this is a more accurate 
measurement of the stable band than that conducted in 
HQ, as we are always forcing with the same values of a 
and v. 

Mass transfer model. The model (|l|) was originally 
presented in Ref. Q, together with a more refined ver- 
sion. A stability analysis of the homogeneous state Xj = 
A shows that ripples are stable (unstable) if /'(A) < 
(> 0). For a convex interaction function, the lower stabil- 
ity boundary X m i n therefore lies at the maximum X marg 
where f'(X marg ) = 0. To investigate nonlinear pattern 
evolution within this model, Eq. (Q) is supplemented by 
the rule that ripples which reach zero length are annihi- 
lated and removed from the system of equations, while 
the remaining ones are relabeled. 

Simulations of the model starting from typical initial 
conditions []l3f in the unstable band A < X marg show that 
an equilibrium wavelength is reached which is essentially 
independent of the initial wavelength, and depends only 
on the interaction function /(A). To gain some insight 
into the selection mechanism, it is useful to recast the 
model into potential form by writing it in terms of the 
position of the troughs between the ripples Xj defined by 



A; 



Then Xj — —dV/dxj, with the potential 



V given by 



A' 



dXf(X), 



(2) 



where N is the number of ripples. It is then plausible to 
conjecture that the equilibrium length A can be found by 
minimizing V for homogeneous ripples, under the con- 
straint that the total length L = NX is conserved. This 
implies that A is determined through the Maxwell con- 
struction applied to /: 



dX /(A) = A/(A). 



(3) 



Comparison with numerical simulations shows that (|^) 
systematically overestimates A, with a better perfor- 
mance the steeper the stable branch of /. Our interpre- 
tation is that the deterministic dynamics gets stuck in 
the multitude of metastable states of (|l|); recall that any 
homogeneous state with A > X marg is a stable, station- 
ary solution. Since the mean ripple length can increase 
only by annihilations, its evolution freezes once all rip- 
ples are in the stable band. The wavelength predicted 
by (|^) should therefore be an upper bound on the actual 
equilibrium wavelength, which is true in all cases we have 
considered. 

Data analysis. Given the time series Xj(t) we want to 
(i) evaluate how well the model ([!]) describes the evolu- 
tion of the ripples and (ii) extract the interaction function 
/(A). For the analysis, we write (Q) in the more general 
form 



Xj = -MXj.,) + 2f c (Xj) - f r (X j+1 ) 



(4) 



with an additional degree of freedom as the functions are 
not required to be equal. 

We want to determine the optimal transformations 
fl(Xj-i), fc(Xj), f r (Xj+i) in the sense that they mini- 



the error X 2 = St 



A, 



^(t)+fl(t)-2f c (t)+f r (t) 

The problem is solved numerically by the ACE algorithm 
EyTsj , The algorithm works by varying fi tC , r m the 
space of all measurable functions until convergence to the 
absolute minimum of \ 2 p~4) . The results are nonpara- 
metric functions which are given in numerical form by 
points, e.g., (Xj, f c (Xj)). An upper bound on the error, 
e.g., on f c , is given by a fc = maxx(df c /dX) ■ <7\, with u\ 
being the measurement error at the point A; analogous 
estimates apply to f fi, f r . The error in the points on 
the A-axis is estimated using the above result and equals 
roughly the errors in the given values of A. 

The quality of each of the resulting functions is given 
by the maximum correlation ipi tC ,r of one of the terms 
with the sum of all the others. A value of ip = 1 implies 
a perfect result, lower values indicate either imperfect 
modeling or (measurement) noise or both. 
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Results. We have performed the analysis on data 
from 18 realizations of a coarsening process initiated with 
approximately 60 small ripples and run with a = 6 cm 
and v — 0.6 Hz. Fig. [I] shows the three functions 
f c and f r averaged over the 18 runs. Clearly there is 
some noise, but the functions are very similar as expected 
from ([l]). The maximum correlations are found to be 
*i = 0.66, * c = 0.88, * r = 0.72. The uncertainty in 
estimating A amounts to about 2 mm, with a maximum 
slope of / of around 3 x 10~ 3 (period) -1 , the absolute 
error of / is 6 x 10 -4 cm/period. This means that the 
result does not fluctuate due to lack of data, rather the 
model (^) does not account for all the variation in the 
ripple evolution. We will return to this point later. 

The maximum of the function lies around X ma rg = 
1.08a, which is a little larger than what was found by 
the numerical simulations in 0], Concerning the shape 
of the function, it is interesting that the slope in the un- 
stable band A < X marg is much larger than the slope in 
the stable band A > X marg . As the slope is a measure of 
the time scale of the dynamics, this implies that the ini- 
tial coarsening stage is much faster than the equilibrating 
stage. 

To use the result function to integrate (|lj) numerically 
for evaluating (^), information about the smallest ripples 
is needed. In lack of such information, we have extrap- 
olated linearly to A = 0, Fig. [I], dashed line. Varying 
the slope of the extrapolated part by a factor of two in 
either direction does not change the final number of rip- 
ples. Fig. |c shows a space-time plot of the numerical 
integration. The model predicts a qualitatively similar 
behaviour as the experiment, however the instability of 
the small ripples develops slower in the model than in 
the experiment. If the extrapolated slope for the interac- 
tion function is made steeper, this instability will evolve 
faster in the model. The model also overestimates the 
final ripple length. For similar initial conditions and sys- 
tem size, the final no. of ripples is typically 18, whereas 
the experiment yields 20. 

Evaluating the Maxwell construction (||) using the 
measured function and extrapolating to account for the 
larger ripples, gives A = 2.2a, which lies far outside of the 
range of definition of /(A). Evidently, the experimentally 
determined /(A) belongs to the class of interaction func- 
tions for which the analytic bound is not useful. 

The upper bound of the stable band, where new rip- 
ples are created, was found to be X max = 1.45a. At 
this point, an innnitesimally small ripple, inserted in the 
trough between two larger ripples, is able to gain mass 
from the neighbours, and thus grows. In the model, this 
corresponds to ripple lengths for which the interaction 
function is smaller than the value of the interaction func- 
tion for ripples at zero length, i.e., f(X max ) = /(0) M. 
Returning to the measured interaction function in Fig.]lj, 
we see that the interaction function at X ma x is approx- 
imately zero, and at A = it is extrapolated to —0.015 
cm/period. Even with the large uncertainty inherent in 
the extrapolation, this does not fit with the theoretical 



picture given above. We can therefore conclude, that the 
model in the form ([l]), is not able to quantitatively pre- 
dict creation of new ripples. The reason for this apparent 
failure relates to the model being developed essentially 
as an expansion around a homogeneous state ||. But 
in the case of a creation, the state is as inhomogeneous 
as it can be: a very tiny ripple flanked by large ripples. 
To account for this, the interaction function should be 
described as a function of the length of the ripple creat- 
ing the separation bubble, Aj, but also of the length of 
its neighbours. Ongoing numerical work indicates that 
it is most relevant to write / = /(Aj, Ai + i)wherc Aj is 
the length of the ripples creating the separation bubble, 
and Aj+i is the one "touched" by the separation bubble 
(see inset Fig. [l]). In fact, the whole coarsening process is 
dominated by highly inhomogeneous configurations, and 
this might also be why simulating ([!]), using the mea- 
sured interaction function, does not produce exactly the 
correct final ripple length. 

Conclusions. We have demonstrated that we can ex- 
tract the interaction function /(A) from spatio-temporal 
data of the evolution of the profile of the sand sur- 
face. The nonlinear data analysis shows its full strength, 
producing nonparametric function estimates, where any 
parametric approach would have failed. With the inter- 
action function it is possible to model the evolution of the 
ripples, in particular their coarsening and equilibration. 

The prediction from the model of the existence of a 
stable band is indeed observed in the experiment. The 
lower bound is predicted at Xmarg /a = 1.08 ± 0.03, while 
the measurement gives X m in/o = 1.13±0.03. In a similar 
experiment || , the lower bound of the stable band was 
found to coincide with the final ripple length: A m ; n = A, 
whereas we find that X m i n < A. The difference between 
the two results might only be apparent, as the number 
of ripples in the annulus of || was 10 or smaller, giving 
rise to an uncertainty in the ripple length on the order of 
the difference between X m i n and A. 

The dynamics responsible for the evolution of a Id rip- 
ple pattern, are also relevant for 2d ripple patterns. How- 
ever, in the latter case, topological defects will be present 
JtJ , and the final length selection might be determined by 
the motion of these || . 
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FIG. 1. The interaction function for a = 6 cm and v — 0.6 
Hz. The thick gray line is a smooth curve drawn through the 
points representing the three different functions. The right 
limit of the plot at A/a = 1.45 is the limit of stability, where 
new ripples are created (X ma x)- The inset shows a sketch 
of the ripples in the part of the oscillation when the flow is 
from the left to the right. The interaction function can be 
interpreted as the transport of sand in the trough, across the 
dashed line. 




FIG. 2. A sketch of the experimental setup seen from above 
(left) and from the side (right). The length of the arm E and 
the width of the channel are not to scale. 
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FIG. 3. a) An example of an extracted profile (dashed line) 
and the fitted triangles with constant slope. The line is shown 
above the profile for clarity, b) The experimental evolution 
of the position of the ripple crests starting from ripples with 
lengths 2.5 cm and evolving with a — 6 cm and ^ = 0.6 Hz. c) 
A simulation of the model (hi) using the extracted interaction 
function and the same initial conditions as above. 
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